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INTRODUCTION 


In  this  report,  we  briefly  describe  the  method  used  and  numerical  results  obtained  for 
diffraction  from  an  absorbing  gold  metallic  grating  with  a  sinusoidal  profile  having  period 

=  2X,  and  =  0.4X,  where  h  is  the  peak-to-valley  ^ting  height  and  X  is  the  wavelength. 
This  grating  has  been  previously  discussed  by  Depine  (Reference  1),  ^d  we  compare  our 
results  below  with  prior  results  presented  in  that  paper.  We  have  previously  compared  our 
results  with  Depine,  as  described  in  Elson  and  Tran  (Reference  2),  where  an  approach  that 
has  similarities  to  that  presented  here  was  used.  Further,  in  this  work  we  will  also  compare 
the  Depine  results  with  yet  another  method  previously  described  by  Elson  and  Tran 
(Reference  3)  where  this  approach  was  also  used  to  investigate  surface  wave  dispersion  on 
truncated  photonic  crystals  (Reference  4).  Finally,  the  method  used  in  this  work  is  identical 
to  that  given  by  Elson  and  Tran  (Reference  5),  which  was  previously  us^  to  cdculate  bulk 
dispersion  of  photonic  media  and  transmission  of  finite  thickness  photonic  media. 

Depine  used  a  surface  impedance  boundary  condition  method.  Elson  and  Tran 
(Reference  2)  used  a  modal  expansion  method  where  Maxwell’s  equations  were  partially 
written  in  Fourier  space  and  the  resulting  differential  equations  were  solv^  as 
eigenfunction  expansions.  The  approaches  shown  in  References  3  and  4  are  very  similar  to 
Reference  2  except  that  Maxwell’s  equations  are  approximated  by  finite  difference 
equations  and  the  solutions  are  obtained  entirely  in  real-space. 

As  also  shown  in  Reference  5,  the  approach  used  in  Ais  work  writes  Maxwell’s 
equations  entirely  in  finite  difference  form  and  then  numerically  integrates  through  the 
profile  region  to  obtain  the  solution.  As  in  References  2  through  5,  the  R-matrix 
propagation  algorithm  is  used  because  of  inherent  numerical  stability.  The  calculation 
method  is  outlined  below. 


CALCULATION  OUTLINE 


We  consider  a  sinusoidal  grating,  shown  in  Figure  1,  which  is  invariant  in  the  y 
direction  and  periodic  in  the  x  direction.  The  profile  shape  is  z  =  fix)  and  is  confined  to  a 
finite  thickness  region  in  the  z  direction.  The  semi-infinite  substrate  and  superstate  are 
homogeneous  and  the  finite  thickness  region  in  between,  which  contains  the  grating 
profile,  is  described  by  a  spatially  variable  permittivity  Eix,z).  The  substrate  (z  <  0)  has 
permittivity  and  the  superstrate  (z >  A)  has  permittivity  (1,0).  When  0<z<h  and 
0  <  z  <  /(a:),  then  e(x,z)  =  When  0<z</iand/i>z>  fix),  then£(x,z)  =  (1,0). 
To  solve  the  diffraction  problem,  we  must  obtain  the  solutions  to  Maxwell’s  equations 
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figure  1.  Schematic  of  Sinusoidal  Grating  Profile  and  Associated  Nomenclature.  The 
finite  difference  increments  are  Ax  and  Az.  A  row  of  =  23  digitized  Ar-points  is  shown 
exphcitly  at  z  =  0.022  micrometers.  There  are  /i^  =  10  divisions  of  height  h  that  yields 

for  this  example  Az  =  0.04A  =  0.022  micrometers  (X,  =  0.55  micrometers).  The  shaded 
regions  within  the  profile  region  show  the  actual  profile  shape  based  on  the  number  of 

layers  n,  since  the  permittivity  £(x,z)  is  independent  of  z  within  each  Az  layer. 


throughout  the  profile  region.  To  do  this,  we  digitize  the  x-coordinate  over  period  d  into 
discrete  x  valuer  where  each  point  is  separated  by  Ax  =  d/n^.  We  also  discretize  the 
z-coordinate  into  values  over  height  h  where  each  point  is  separated  by  Az  =  h/n^. 
Maxwell’s  equations  are  now  approximated  with  finite  difference  expressions: 
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V  xE(a:,z)  =  i(fO/c)B(x,z): 


E^{x,z  +  ii>z)-E^{x,z) 
Az 


+- 


ic _ \By{x  +  liiX,z)-By{x,z)  ^  By{x-^x,z)-By{x,z)\ 


(o{^xy 


e(jc  + Ax/2,z) 


£{x-^xxll,z) 


(la) 


E  {x,z^^)-Eix,z)  io)  .  ^ 

— - - - - - = - B^{x,z  +  Az) 

Az  c 

V  X  =  -i(0)/c)£{x)E{x,z): 


(lb) 


BJx,z  +  Az)  -  B^{x,z)  ico 


Az 


= - £(x,z)E  {x,z) 

c 


+-^^^^^[2Eyix,z)  -  Eyix  +  Ax,z)  -  Ey(x  -  Ax:,z)} 


(Ic) 


Byix,z  +  Az)-By(x,z) 
Az 


=  —  £ix,  z  +  Az)£';t  (x,  z  +  Az) 
c 


(Id) 


In  Equation  1,  x  and  z  are  each  a  discrete  coordinate  in  the  set  of  n,  and  points.  These 
equations,  which  apply  to  the  profile  region,  are  valid  from  z  =  0  to  z  =  /i.  For  numerical 
purposes,  we  truncate  the  problem  so  that  the  integer  index  ranges  from  -N  -»  N  and 
this  means  that  when  all  x-points  are  included.  Equation  1  represents  4(2A^+1)  equations 
for  a  given  discrete  z- value.  If  we  include  all  2N  +  1  x-points.  Equation  1  may  be  arranged 
into  a  matrix  equation  having  the  form 


^  E,iX,z)  " 

^  5,(X,z)  > 

Ey(X,z) 

By{X,z) 

£,(X,z  +  Az) 

~  *2 

B,(X,z  +  Az) 

^E/X,z  +  Az)^ 

^ByiX,z  +  Az)y 
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where  r,  and  Fj  are  square  matrices  obtained  from  Equation  1  and  X  represents  the  set  of 
discrete  jc-points.  The  columns  containing  the  x  and  y-components  of  E  and  B  are 
column  vectors  where  each  of  the  E  and  B  components  have  2A^+1  elements  and  each 
element  denotes  a  different  j:-coordinate.  Defining  a  matrix  r  =  r,~'r2  and  abbreviating  the 
column  vectors  with  tilde  notation  as 


E  = 

k) 

and  B  = 

k) 

yields  Equation  2  as 


'  E(X,z)  ' 

=  r(Az) 

"  B(X,z)  ^ 

^E(X,z  +  Az)^ 

^B(X,z-t-Az); 

(3) 


(4) 


This  equation  defines  the  r  matrix  for  the  layer  bounded  by  z-^z  +  Az  and  needs  to  be 
calculated  for  each  of  the  layers.  In  constructing  the  r  matrix,  certain  terms  of 
Equations  la  and  Ic  require  special  attention.  We  note  that  for  extreme  x-values  within  a 
period,  some  field  terms  with  arguments  x±Ax  will  fall  outside  the  dimension  of  one 
period.  These  terms  are  “wrapped  around”  by  using  periodic  boundary  conditions  (Bloch 
theorem).  As  described  in  detail  in  References  2  through  5,  the  fields  at  the  z  =  0  and 
z  =  li  planes  may  be  related  by  a  global  relationship  as 


^E(X,0)' 

,E(X,/i), 


=  R(ft) 


^B(X,0)^ 

^B(X,/I); 


(5) 


where  the  square  matrix  R(h)  is  obtained  by  a  recursive  algorithm  involving  the  r  matrices 
for  each  of  the  layers.  Equation  5  is  described  in  real-space,  and  it  is  convenient  at  this 
time  to  transform  to  Fourier  space  by  means  of  the  square  matrix  F(K,X)  where 


F(K,X) 


^E(X,0)^ 

^E(X,/t), 


^E(K,0)^ 

,E(K,/z), 


(6) 


and  similarly  for  the  magnetic  field.  The  vector  K  denotes  the  set  of  wave  vectors  in 
Fourier  space  analogous  to  the  set  X  of  x-coordinates  in  real  space.  We  apply  this  F 
matrix  to  Equation  5  that  yields 
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^E(K,0)^ 

=m) 

^B(K,0)^ 

^E(K,;i); 

,B(K,ft), 

(7) 


where  the  matrix  R(/i)  =  FR(ft)F  * .  Since  the  z  —  0  and  z  —  h  planes  define  the 
boundaries  of  the  profile  region  adjoining  the  substrate  and  superstrate,  respectively,  we 
may  match  boundary  conditions  at  these  boundaries  with  the  respective  fields  in  the 
substrate  and  superstrate.  At  the  substrate  boundary,  the  boundary  conditions  may  be 

written  as  E(K,0)  =  E'(K,0)  and  B(K,0)  =  B'(K,0)  where  the  superscript  t  refers  to  the 
transmitted  field  region  (substrate).  At  the  superstrate  boundary,  the  boundary  conditions 
are  E(K,/t)  =  E'”^(K,/i)  +  E"(K,/i)  and  B(K,A)  =  B'""(K,/t)  +  B''(K,/i)  where  the 
superscripts  inc  and  r  refer  to  the  incident  field  and  reflected  field  region  (superstrate), 
respectively.  Using  these  boundary  conditions  in  Equation  7  yields 


E'(K,0) 
E"’"(K,/i)  +  E7K,/i)y 


\  _  f 

=m) 


B'(K,0) 
B'^(K,/i)  +  r(K,/i) 


(8) 


While  transformation  to  Fourier  space  is  convenient  for  diffraction  calculations,  it  is  also 
important  since,  in  the  homogeneous  superstrate  and  substrate  regions,  a  relation  between 
E  and  B  can  be  written  as 


r(K,z)  =  Z^'E'(K,z) 


where  j  —  roTt.  The  form  of  the  Z  is  given  in  References  2  through  4.  With  this 
relationship,  we  may  finally  write  the  end  result  as 


-R„Z' 


l-RjjZ^ 


Y  E'(K,0)^ 

,W(K,h)j 


^0 

.-I 


R,2  Y  E'""(K,/i)^ 
R22AB'"^(K,/i) 


(9) 


where  the  square  matrix  R  has  been  split  into  four  quadrants  each  denoted  by  R;„„  with 
the  appropriate  m  and  n  values.  This  matrix  equation  may  be  solved  by  standard  means  to 
calculate  the  transmitted  E'  and  reflected  E'^  fields. 
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NUMERICAL  RESULTS 


We  assume  a  plane  wave  of  wavelength  X  =  0.55  micrometers  is  incident  on  the 
sinusoidal  grating  shown  in  Figure  1,  where  the  incident  polarization  may  be  perpendicular 
or  parallel  to  the  grooves.  The  grating  profile  is  given  by  /(jc)  =  /i[l  +  cos(2;a:/(i)]/2. 

The  gold  substrate  material  has  permittivity  e^„i,  =  (-5.28,  1.48)  atX  =  0.55  micrometers. 
The  numerical  results  are  given  for  polarization  parallel  to  the  plane  of  incidence 
(perpendicular  to  the  grooves).  The  numerical  parameters  used  include  =  83  jc-points 
which  is  the  same  as  the  number  of  diffracted  orders  and  this  yields  Ax  =  0.0241 A  = 
0.0136  micrometers.  We  let  Az  =  O.lAx  which  yields  =  165. 

In  Figure  2  we  compare  the  present  calculation  method  with  the  methods  described 
in  References  2  through  4.  The  input  parameters  and  for  the  respective  method  are 


FIGURE  2.  Comparison  of  Difftaction  Intensity  Versus  Angle  of 
Incidence  For  Three  Methods  of  Calculation.  The  solid  line  is  for  the 
present  numerical  integration  method.  The  dotted  line  is  for  the  yt-space 
method  described  in  Reference  2.  The  dash  line  is  from  the  real-space 
method  described  in  References  3  and  4.  The  curves  to  the  left  and  right 
of  the  vertical  dotted  line  are  for  the  +1  and  -3  order,  respectively.  It  is 
seen  that  the  three  methods  are  in  close  agreement. 
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shown  in  References  2  through  4,  but  the  grating  design  and  material  parameters  are  the 
same  as  used  here.  It  is  seen  that  the  three  methods  agree  quite  well. 

Figure  3  compares  the  Depine  Reference  1)  data  with  the  present  numerical  integration 
method.  There  is  considerable  disagreement  and  the  reason  is  not  clear.  In  Depine’s 
paper,  he  shows  some  numerical  data  that  is  based  on  a  rigorous  integral  equation  method 
and  these  data  agree  quite  weU  with  his  approximate  impedance  boundary  condition 
approach. 


FIGURE  3.  Diffraction  Intensity  Versus  Angle  of  Incidence.  These  data 
compare  the  Depine  results  (dotted  line  with  solid  circles)  with  the  present 
numerical  integration  method  (solid).  The  vertical  dash  line  separates  the  +1 
and  -3  diffracted  orders.  It  is  seen  that  there  is  moderate  disagreement. 


CONCLUSIONS 


We  have  described  a  method  of  calculating  grating  diffraction  that  has  the  advantage  of 
not  requiring  calculation  of  eigenvectors  and  eigenvalues.  A  disadvantage  is  that  the 

numerical  integration  through  the  profile  (z-direction)  requires  that  the  increment  Az  «  X,. 
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The  method  is  versatile  in  that  different  grating  profiles  are  simple  to  incorporate  and 
extension  to  multilayered  grating  stmctures  is  possible  in  analogy  with  References  2  and  3. 

The  numerical  results  for  the  numerical  example  chosen  here  show  agreement  of  the 
present  method  with  two  prior  methods.  Differences  between  these  three  methods,  which 
are  evident  in  Figure  2,  may  be  decreased  by  varying  discretization  (numerical)  parameters 
which  can  affect  accuracy.  The  present  method,  and  hence  all  three  methods,  are  in 
moderate  disa^eement  with  Depine  (Figure  7  contained  in  Reference  1).  The  reason  for 
disagreement  is  left  to  future  study. 
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